#include <algorithm>
#include "Core/numlib.h"
#include "Embedded/EmLevelOp.h"

template <int Dim>
Real EmLevelOp<Dim>::NRM(const Tensor<Real, Dim> &aData, int staggered, int p) const
{
  if(staggered == CellCentered) {
    const auto dx = emGrid.spacing();
    const auto &nodeLabels = emGrid.nodeLabels;
//    using NodeType = typename EmGrid::NodeType;
    if(Dim >= 3) {
      Real r = 0.0;
      if(p == 0) {
#pragma omp parallel for default(shared) schedule(static) reduction(max:r)
        loop_box_3(emGrid, i, j, k) {
          if(nodeLabels(i, j, k) >= 0)
            r = std::max(r, std::abs(aData(i, j, k)));
        }
        return r;
      } else if(p == 1) {
#pragma omp parallel for default(shared) schedule(static) reduction(+:r)
        loop_box_3(emGrid, i, j, k) {
          if(nodeLabels(i, j, k) >= 0)
            r += std::abs(aData(i, j, k));
        }
        return r * prod(dx);
      } else if(p == 2) {
#pragma omp parallel for default(shared) schedule(static) reduction(+:r)
        loop_box_3(emGrid, i, j, k) {
          if(nodeLabels(i, j, k) >= 0)
            r += ipow<2>(aData(i, j, k));
        }
        return sqrt(r * prod(dx));
      }
    } else if(Dim >= 2) {
      Real r = 0.0;
      if(p == 0) {
#pragma omp parallel for default(shared) schedule(static) reduction(max:r)
        loop_box_2(emGrid, i, j) {
          if(nodeLabels(i, j) >= 0)
            r = std::max(r, std::abs(aData(i, j)));
        }
        return r;
      } else if(p == 1) {
#pragma omp parallel for default(shared) schedule(static) reduction(+:r)
        loop_box_2(emGrid, i, j) {
          if(nodeLabels(i, j) >= 0)
            r += std::abs(aData(i, j));
        }
        return r * prod(dx);
      } else if(p == 2) {
#pragma omp parallel for default(shared) schedule(static) reduction(+:r)
        loop_box_2(emGrid, i, j) {
          if(nodeLabels(i, j) >= 0)
            r += ipow<2>(aData(i, j));
        }
        return sqrt(r * prod(dx));
      }
    }
  } else {
    assert("To be implemented.");
  }
  return 0.0;
}

template <int Dim>
void EmLevelOp<Dim>::copyGhosts(const Tensor<Real, Dim> &aSrc, Tensor<Real, Dim> &aDest) const
{
  for(const auto &idx : emGrid.ghostNodes)
    aDest(idx) = aSrc(idx);
}

template <int Dim>
void EmLevelOp<Dim>::zeroExtNodes(Tensor<Real, Dim> &aData, bool zeroGhosts) const
{
  const auto &nodeLabels = emGrid.nodeLabels;
//  using NodeType = typename EmGrid::NodeType;
  ddfor(emGrid, [&](const iVec &idx) {
    if(nodeLabels(idx) < 0)
      aData(idx) = 0.0;
  });
}

template <int Dim>
void EmLevelOp<Dim>::zeroExtNodes(Tensor<Real, Dim + 1> &aData, bool zeroGhosts) const
{
  auto comps = getComp(aData.box(), Dim);
  for(int k = comps.lo()[0]; k <= comps.hi()[0]; ++k) {
    Tensor<Real, Dim> alias(aData, k);
    zeroExtNodes(alias, zeroGhosts);
  }
}

//============================================================
template class EmLevelOp<SpaceDim>;
